# 龙贝格积分

import numpy as np


def romberg(f, a, b, n):
    r = np.zeros((n, n), dtype=float)  # romberg 下三角矩阵
    ba = b - a  # 区间总长度
    h = np.array([(ba / 2**i) for i in range(0, n)])  # 步长列表
    r[0, 0] = ba * (f(a) + f(b)) / 2.0

    for j in range(1, n):
        # 计算新的 第一列
        subtotal = 0.0
        for i in range(0, 2**(j - 1)):
            subtotal += f(a + (2 * i + 1) * h[j])
        r[j, 0] = 0.5 * r[j - 1, 0] + h[j] * subtotal
        #  计算新的 2,...,j 列
        for k in range(1, j + 1):
            r[j, k] = (4**k * r[j, k - 1] - r[j - 1, k - 1]) / (4**k - 1)
    return r


def test():
    np.set_printoptions(formatter={'float': '{: =-#.14f}'.format})
    f = lambda x: np.log(x)
    r = romberg(f, 1, 2, 4)
    print(f'romberg matrix is:\n{r}')


if __name__ == '__main__':
    test()
